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Abstract  - The  accurate  and  fast  evaluation  of  surface  reaction  integrals  for  Method  of  Moments  computations  is 
presented.  Starting  at  the  classification  of  the  integrals  into  regular  and  weakly,  strongly  and  nearly  singular 
integrals,  appropriate  methods  are  presented  that  handle  each.  A Gauss-Legendre  quadrature  rule  evaluates  regular 
integrals.  For  singular  integrals,  the  singularity  is  lifted  or  weakened  by  an  extraction  of  the  singularity,  a transform 
to  polar  coordinates  or  a domain  transform.  The  resulting  regular  integral  is  in  turn  solved  by  a quadrature  rule.  The 
different  methods  are  finally  applied  to  an  example,  and  the  resulting  accuracy  tested  against  the  analytical  result.  The 
presented  methods  are  general  enough  to  be  used  as  integration  methods  for  integrands  with  various  degrees  of 
singularity  and  is  not  limited  to  Method  of  Moments. 

1.  INTRODUCTION  comparison  of  the  methods  closes  the  presentation. 


The  solution  of  integro-differential  equations  by  the  Method 
of  Moments  is  directly  influenced  by  the  accuracy  of  the 
matrix  elements  corresponding  to  the  reaction  integrals. 
Efficient  and  accurate  methods  for  the  evaluation  of  the 
reaction  integrals  are  hence  necessary.  The  numerical 
methods  must  fulfill  two  major  requirements.  First,  the 
method  should  converge  quickly  to  ensure  a low 
computation  time.  Second,  the  evaluation  should  be  accurate 
in  order  to  obtain  a good  overall  accuracy. 

The  evaluation  of  the  reaction  integrals  gives  rise  to 
integrals  whose  kernel  include  the  free-space  Green’s 
function  and  its  derivatives  [PM73].  This  necessitates  the 
proper  treating  of  four  kinds  of  integrands.  First,  regular 
integrals  are  solved  by  usual  quadrature  rules.  Second, 
weakly  singular  integrands  require  mathematical  transforms 
before  integrating.  Third,  strongly  singular  integrands  are 
only  accesible  to  specialized  integration  methods.  Finally, 
observation  points  that  are  located  very  close  to  the 
integrational  domain  give  rise  to  nearly  singular  integrals. 
This  paper  presents  powerful  integration  methods  for  each 
integral. 

The  first  purpose  of  the  paper  is  to  present  different  methods 
and  their  applicability  to  the  different  types  of  integrals.  The 
second  purpose  is  to  present  a new  accurate  method  to 
evaluate  strongly  singular  kernels.  The  third  purpose  is  to 
carefully  compare  the  appropriate  integration  methods.  To 
the  authors'  knowledge,  the  present  complete  description  of 
the  integration  methods  is  not  available  in  the  literature  as 
such. 

The  organization  of  the  paper  is  as  follows.  The  following 
section  presents  the  classification  of  the  kernels.  Then,  the 
four  methods  for  regular,  nearly  singular,  weakly  singular 
and  strongly  singular  integrands  are  presented  with  an 
example  of  a rectangular  surface  element  [HvHWOO].  A 


2.  KERNEL  CLASSIFICATION 


The  reaction  integrals  within  a Method  of  Moments 
computation  involve  the  integration  of  the  Green's  functions 
or  their  derivative  together  with  the  basis  functions  and 
testing  functions  [PM73]. 

The  free-space  Green's  function  with  a source  point  at  r' 
G(F,r')  = e'mun/{An\r  -r'§  and  its  derivative  VG(F,F') 
= -V'G(F  ,F ')  determine  the  behavior  of  the  kernels.  The 
distance  factor  R contained  in  the  Green's  function  is: 

r = r - r\  = - *7 + (/  - yf + (*'  - z')2  a) 


and  determines  if  the  integrand  at  the  point  P is  regular  (for 
r * r ' ) or  singular  (for  r =*  r ' ). 


For  the  following,  the  propagation  constant  ft  is  separated 
into  a real  and  an  imaginary  part: 

P = o)^fi0(s'-je"-jK/(o)  =P'  + jP"  (2) 


with  the  usual  material  parameters  e , ^ and  k and  the 
angular  frequency  w . 

One  obtains  then: 


-jpR 


P"R 


4 nR  4:zR 


(cos(/3'/?)-ysin(^'/f)) 


(3) 


The  derivative  V'G(r,r ')  in  cartesian  coordinates  is  then: 


V'G : 


- / - p"R 

r - r e 
R3  4jt 


j^-cos(/3'/?j  + ysin(/37?) 


+ 


P"R 


r -r  e 

R*  An 


[-/}'  sin(/37?)  + P"  cos(/37?) 
-j(p'cos(p'R)  + P"sm(p'R))] 


(4) 
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Consider  each  term  in  (3)  and  (4)  when  the  observation 
point  P moves  closer  to  the  source  point,  or  R -*•  0 . In  this 
case,  (3)  and  (4)  become  singular,  the  order  of  their 
singularity  is  listed  in  Table  1.  The  imaginary  part  (5)  of 
the  Green’s  function  is  regular,  and  further,  the  derivative  of 
the  Green's  function  has  a regular  term  (6).  Mathematically, 
(5)  and  (6)  can  be  integrated  with  no  further  problems. 

The  real  part  of  the  Green's  function  (7)  and  parts  (8)  of  the 
derivative  of  the  Green's  function  are  weakly  singular  with 
terms  - 1 /R  . This  singularity  can  be  easily  regularized. 

The  last  term  (9)  caused  by  the  derivation  of  the  Green's 
function  is  strongly  singular  ~1//?2  for  R-*  0.  A 
regularisation  must  be  carried  out  before  an  integration  can 
be  carried  out.  For  the  strongly  singular  case,  the 
regularisation  is  more  complicated  than  for  the  weakly 
singular  case.  Appropriate  integration  methods  are 
necessary. 

3.  Integration  Methods 

The  integration  methods  for  the  kernels  in  Table  1 are 
chosen  according  to  the  following  three  cases: 

- R > d:  the  distance  between  the  observation  point  P and 
the  source  point  is  greater  than  a predetermined  minimal 
distance  d (regular), 

- R < d:  the  distance  between  the  observation  point  P and 
the  source  point  is  smaller  than  a predetermined  minimal 
distance  d (regular,  but  nearly  singular), 

R — * 0 : observation  point  and  source  point  P coincide 
(weakly  and  strongly  singular). 

Table  2 shows  a summary  of  which  integration  method  is 
applicable  for  the  various  source-observation  point 
distances.  As  seen,  special  integration  methods  are  needed 
that  correctly  treat  regular,  nearly  singular,  weakly  singular 
and  strongly  singular  integral  kernels.  They  are  in  detail: 

- Regular  integration : Solved  by  a Gauss-Legendre 
quadrature  method.  The  method  is  numerically  accurate 
and  easy  to  implement.  Section  3.1  sketches  the  method. 

- Nearly  singular  integration:  For  a distance  R smaller 
than  a minimal  distance  d , the  integrand  changes  rapidly 
and  becomes  nearly  singular.  The  appropriate  integration 
method  is  laid  out  in  section  3.4. 

- Weakly  singular  integration:  The  method  for  weakly 
singular  kernels  is  outlined  in  section  3.2.  The 
integration  involves  either  a domain  transform,  a 
transform  to  polar  coordinates  or  an  extraction  of  the 
singularity. 


Table  1 Classification  of  the  integrands  for  R 0 


regular  kernels 

- <5> 

‘Me,'R)- ,fs- -MfrR)]  (6) 


weakly  singular  kernels 


- UR 


*r* 


^—cos(p ,R) 
4jiR 


r’-re 


r 


R2  4 Jt 


P”  cos  {p'*)-j\p  cos  ( P'R )- 


sin(/3  ’R) 


R 


(7) 


(8) 


strongly  singular  kernels 


f'-r  S'* 
Ry  4/r 


cos 


(,37?) 


(9) 


Table  2 Integration  methods  for  the  various  source- 
observation  point  distances  R 


regular  kernels 

R>d) 

R <d  > regular  integration 

-►  section  3.1 

R-*  oJ 

weakly  singular  kernels  ^ 

R > d : regular  integration 

section  3.1 

R < d : nearly  singular  integration 

section  3 .4 

R 0 * weakly  singular  integration 

— ► section  3 .2 

singular  kernels 

R > d : regular  integration 

-»  section  3.1 

R < d : nearly  singular  integration 

-*  section  3 .4 

R — ► 0 : strongly  singular  integration 

section  3 .3 

- Strongly  singular  integration:  Two  methods  in  section 
3.3  carry  out  the  integration  of  strongly  singular 
kernels.  The  methods  are  composed  of  the  methods  for 
weakly  singular  and  regular  kernels.  Section  3.3.2 
presents  a new  method  that  combines  the  domain 
transform  and  the  extraction  of  singularity.  Transforms 
and  an  extraction  of  singularity  regularize  the 
singularities,  the  resulting  kernels  are  integrated  by  the 
Gauss-Legendre  quadrature  as  in  3.1 . 
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3.1  Regular  Integration 


A direct  Gauss-Legendre  quadrature  is  proposed  for  the 
regular  parts  in  Table  2 [Abr70],  [Sch93j.  In  (10),  the 
general  Gauss-Legendre  quadrature  for  a general,  regular 
integrand  f(w,,w2)  is  given  as: 


JJ Fiu>  >«2  ){g{u,’ui  )dutdu2 


(10) 


F is  a function  of  two  continuous  variables  ux  and  u2.  They 
correspond  to  the  surface  variables  on  a surface  element  in 
the  Method  of  Moments  [HvHWOO].  The  expression 

^8{ui>ui)du\du2  is  the  differential  surface  element  that 

expresses  the  transform  of  the  surface  differential  into  the 
Mj,m2  system. 


For  the  numerical  integration,  the  « , u}  are  the  nodes  on 
the  surface  in  the  interval  =[-l,+l].  The  w.  and  w.  are 

the  weights  for  the  nodes.  Both  are  listed  in  the 
mathematical  literature  [Abr70]  and  are  not  repeated  here. 


3.2.1  Transform  to  Polar  Coordinates 

The  singular  point  P is  arbitrarily  located  on  the  surface 
element,  with  its  normalized  coordinates  ^*[-1,1], 
u2  - [-1J]  • The  weakly  singular  integrands  are  regularized 
by  a transform  into  polar  coordinates  [GG90].  The  process 
of  regularization  is  depicted  in  Fig.  2.  First,  the  surface 
element  is  partitioned  into  four  rectangles  Vk  with  the 

singular  point  P located  on  the  comer  that  is  common  to  all 
rectangles.  Each  rectangle  is  then  in  turn  subdivided  into 
two  triangles,  the  integration  is  carried  out  on  each  of  the 
eight  triangles  separately. 


Figure  2 Sub-division  of  a discretization  element  into 
four  quadrangles  Vk  and  subsequent  transform 
to  polar  coordinates. 


3.2  Weakly  Singular  Integration 

Three  different  methods  are  proposed  in  this  section  for  the 
numerical  computation  of  weakly  singular  kernels  if 
0: 

- section  3.2.1:  regularization  of  the  kernel  by  a transform 
into  polar  coordinates 

- section  3,2.2:  regularization  by  a domain  transform 

- section  3.2.3:  regularization  by  an  extraction  of  the 
singularity  and  subsequent  analytical  integration 

All  three  integration  methods  are  presented  for  the  sake  of 
completeness.  The  methods  are  later  also  used  for  the 
integration  of  strongly  singular  kernels. 

As  an  example  consider  the  weakly  singular  term 
(x'-xf)/{4xR2)P'  cos (P'R)  (8)  on  a planar  rectangular 
element  as  in  Fig.  1 with  a size  of  X/10  which  is  typical  for 
a surface  element  in  a Method  of  Moments  implementation. 


The  transform  into  polar  coordinates  p , ft  is 
Wj  +pcostf  (11) 

ui  m ui  + psintf  (12) 

duxdu2  « pdpdft  (13) 

The  integral  is  then  in  polar  coordinates 

+i  +i 

f f F(u,  >“2  )^ut,u2)duidu1 

(14) 

■sr*  ( #.p(  # ( #,p(  i jjjpdpdft 

The  weight  p that  arises  due  to  the  transform  from  the 
uiyu2 -system  into  the  p,#-system  lifts  the  singularity  and 
regularizes  the  integrand.  The  numerical  integration  of  (14) 
is  carried  out  by  observing  that  the  upper  limit  p/[ft)  is  a 
function  of  ft. 


■*P=  “T  = 0.25m 


Figure  1 Example  geometry  to  outline  integration 
accuracy. 


For  the  above  example,  the  regularization  of  the  integrand  is 
shown  in  Fig.  3.  Fig  3a)  shows  the  weakly  singular 
integrand  in  the  ux  ,«2 -system.  Fig  3b)  to  3e)  depict  the 

regularized  integrand  in  polar  coordinates.  The  angle  ft  is 
defined  on  the  following  intervals  for  the  triangles  in  3.2.2: 

- Fig.  3b):  4,02  s # < 5,18  and  /*{#)  - -l,5/sin(#), 

- Fig.  3c):  -1,1  s ft  < 0,59  and  p/[d)  - 0,75/cos(tf), 

- Fig.  3d):  0,59  s # < 2,76  and  /y( fl)  - 0,5/sin(tf) , 

- Fig.  3e):  2,76  s < 4,02  and  /&(#)  - -l,25/cos(#). 
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duxdu2  « 2 ADtdsxds2 


4 


dtxdt2 


(19) 


ADl  are  the  surfaces  of  the  triangles  in  the  w,  ,w2 -system. 
For  source  and  observation  points  that  are  not  in  the  centre 
= 0 , u2  - 0 of  a surface  element,  the  ADi  are  different  for 
each  triangle. 


Figure  3 Transform  of  a weakly  singular  integrand  into 
regular  sub-integrands  using  a transform  into 
polar  coordinates.  Each  integrand  is  then 
normalized. 


As  seen  in  Fig.  3 b)  to  e),  the  integrands  present  a smooth 
behavior. 


3.2.2  Domain  Transform 

Alternatively  to  the  previous  transform,  the  following 
domain  transform  may  be  used.  First,  the  rectangle  in  the 
Wj,m2  -system  is  subdivided  (cf.  Fig.  4)  into  the  triangles 
Dj , D2 , D^  and  D4  [Duf82,  Dom93].  The  singular  point  P 
is  then  located  on  the  comer  that  is  common  to  all  triangles. 


Figure  4 Subdivision  of  a surface  element  into  triangles, 
and  subsequent  transform  of  a sub-triangle  Dk 
into  a rectangle. 

A non-linear  domain  transform  maps  each  triangle  into  a 
rectangle,  the  singular  point  P is  mapped  to  a singular  edge. 
The  original  variables  «,  ,w2  are  transformed  to  new  surface 


variables  s{  «[0,l]  and  s2  = [0,l]  by: 

- (l  - •*,  K"  +*,(!-  *2  K‘  ‘ + VA°"  (15) 

u2  = (l  - ^ + sx  (l  - s2  )u2tx  + (16) 

The  subscripts  in  the  DkJ  correspond  to  the  number  of  each 

triangle  k = 1 ....4  and  its  comer  i = 1 ....3.  Now,  the  interval 
is  shifted  to  obtain  symmetric  integration  limits 

*,-(*,+ 1)/2  (17) 

*2»('2+l)/2  (18) 


tx  and  t2  are  now  from  -1  to  +1.  The  differentials  are 
transformed  by  the  Jacobian: 


A weakly  singular  integral  is  hence  regularized  using  the 
integration  in  (20): 


+i  +i  

,«2  )Jg(ul,u2)duldu1 


4 


(fi  +l)^° 


dttdt2 


(20) 


The  7^0  are  the  weights  arising  by  the  transform  of 
the  surface  differential  into  the  «,  ,«2 -system.  The  transform 
of  the  differentials  from  the  w,  ,m2 -system  into  f,  >t2 -system 
is  according  to  (19)  adds  the  term  (*,  +l)  that  lifts  the 
singularity  and  regularizes  the  integrand. 


In  Fig.  5,  the  triangles  arising  in  the  regularization  are 
depicted.  5b)  to  5e)  show  the  regularized  integrands  in  the 
tl  d2 -system.  Again  the  integrands  are  very  smooth. 


Figure  5 Transform  of  a weakly  singular  integrand  into 
regular  sub-integrands  using  a domain 
transform  The  sub-integrands  are  normalized. 


3.2.3  Singularity  Extraction 

The  singularity  of  the  Green's  function  can  be  lifted  by 
subtracting  an  analytically  integrable  part.  The  method 
extracts  the  singularity  in  such  a way  that  the  extracted  term 
is  analytically  integrated  as  in  (21)  for  a singular  point  P: 


By  subtracting  the  inverse  distance,  the  first  integrand  in 
(21)  becomes  regular.  This  term  is  integrated  using  the 
usual  Gauss-Legendre  integration  as  in  section  3.1.  The 
remaining  singular  integral  in  l/R  is  integrated  analytically. 
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For  a linear  triangular  surface  patch,  the  integration  of  HR 
can  be  carried  out  analytically  with  no  approximation.  For 
higher  order  geometrical  surface  approximation  functions 
such  as  the  biquadratic  [HvHWOO],  however,  the  integration 
cannot  be  carried  out  analytically.  In  this  case  the  integrand 
is  developped  into  a series.  For  the  ^-coordinate,  one 
obtains: 


3.3.1  Transform  to  Polar  Coordinates  and 
Extraction  of  Singularity 

The  transform  to  polar  coordinates,  origin  at  the  point  P, 
analogous  to  section  3.2.1,  (11)  and  (12),  ux  - + p cos  , 

and  u2  = u[  + psin  reduces  the  order  of  the  singularity  in 

(9)  by  one  by  virtue  of  the  integration  weight 
duxdu2  * pdpdd' . The  resulting  singularity  is  only  weak, 

and  the  method  for  the  extraction  of  the  singularity  in 
section  3.2.3  is  applicable.  The  term  that  determines  the 
singularity  in  (9)  is  given  by: 


y and  z are  analogous.  Completing  and  combining  the 
expressions  in  xyy , and  z leads  to  : 


R 


■~1  / 


-“if  + 

g22 («,  -«;)'  + 


(23) 


(x1 -xp)/(4nR3)  (25) 

For  an  analytical  integration  in  p,  the  integrand  is  developed 
into  a series  around  the  point  P.  With  (11)  and  (12), 
ux  = p cos  , u2  - u2  - p sin  & , one  obtains  for  an  x- 
component  in  (22): 


-<)  (x'-xf)=p 

dxf 

„ dx' 

COS  17  + 

sini? 

du{ 

Mi-«r 

, in  (23)  is  carried  out 

“»-«T 

u2mU* 

-0(p>) 


(26) 


analytically  [SC95].  With  the  Taylor  series,  one  finally 
obtains: 


sH 


duAu , 


= pA,(#)  + C>(p2) 


Analogous  to  (26),  the  expressions  for  the  y-  and  z- 
components  yield  expressions  in  and  With 

these,  one  obtains  the  term  that  is  necessary  to  express  the 
distance  R: 


(24) 


A(  fl)  = (27) 

The  distance  R depends  now  on  polar  coordinates  and  is 
given  with  (26)  and  (27): 

/?2(p,tf)  = p2AJ(#)  + 0(p3)  (28) 

With  (28)  and  (*’  - xp)/R  = A,(tf)/A(tf)  + O(p)  one 
obtains: 


3*3  Strongly  Singular  Integration 

For  the  integration  of  strongly  singular  integrands,  the 
methods  of  section  3.2  are  combined.  The  combinations 
yield  regular  integrands  that  are  dealt  with  the  method 
described  in  3.1. 

Two  methods  lift  the  singularity: 

- The  application  of  a transform  to  polar  coordinates 
[PG89],  [GG90],  [GC87]  and  subsequent  extraction  of 
the  singularity.  This  is  presented  in  section  3.3.1. 

- A new  method  in  section  3.3.2  combines  the  domain 
transform  and  a subsequent  extraction  of  the  singularity. 
This  new  method  yields  very  smooth  integrands  that  are 
accessible  to  numerical  integrations. 

As  an  example,  consider  the  integration  of  the  singular 
kernel  (9)  in  the  geometry  in  Fig.  1 . 


R}  { A{&)  [P>)  p2A2(i^)  + 0(p3) 

1(AW,  ojp)  ) 

P2  l A3(tf)  A2($)  + 0(p)J 


Incorporating  the  integration  weight  p for  the  transform  of 
the  differentials  from  the  Wj,«2-system  into  the  p,tf- 
system,  one  obtains  the  expression: 


4 JtR* 


duxdu2 


1 

4^rp 


(30) 


The  extraction  of  the  singularity  in  (30)  according  to  (9) 
produces  the  integrand: 


i AW 

4 n pA3( &) 


(*'■ 


■xr)eT 

R3  An 


pcos(/J'r) 


1 AW 

pA3{&) 
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As  the  limit  0[p)jp  O(l)  in  the  first  parenthesis  in  (31) 

exists,  the  singularity  in  the  first  term  is  lifted  and  the 
integrand  is  regular.  The  integration  is  then  carried  out  e.g. 
by  the  method  in  section  3.1  in  polar  coordinates  p( tf)  and 
for  all  rectangles  Vk : 

p A3(#) 

(s'(p,ff)-sf) 


\n 


Jff 

ft.,  o 


R'{p,V) 


’'pcos()37?(p,t7)) 


dpd& 


For  the  evaluation  of  the  extracted  term  in  (31),  the  Cauchy- 
limit  is  computed: 

-j-  lta 


dpdft 


pA3(tf) 


(33) 


dtf 


(P.(#)J 

Now,  introduce  an  e vicinity  using  polar  coordinates: 

-p >M  #)A2(#)  + 0(pt3(#)) 

This  is  then: 

(35)  in  (33),  one  obtains 


(34) 


(35) 


ln(p/„{V)A(V))d& 


with  (see  also  [GC87],  [GG90],  [HRHR97],  [PG89]) 
1 


Atc 


-lim 


(in  e)  = 0 


(36) 

(37) 


The  expression  in  brackets  vanishes  when  carrying  out  the 
integration  and  summation  before  taking  the  limit.  Then  the 
limit  does  not  need  to  be  carried  out  anymore.  The  strongly 
singular  part  in  (9)  caused  by  the  derivative  of  the  Green’s 
function,  is  now  readily  treated  by  (31),  (32),  and  (36). 
Consider  the  integrand  in  (32).  In  Fig.  6a),  the  integrand  is 
shown  in  the  w,  ,w2-system,  It  is  then  transformed  into  the 

sub-integrands  as  in  Fig.  6b)  to  e).  The  sub-integrands  are 
formed  by  the  sub-division  process  in  Fig.  2.  The 
integration  domains  are: 

- Fig.  6b):  4,02  =s  & < 5,18  and  p(tf)  = -l,5/sin(#), 

- Fig.  6c):  -1,1  s < 0,59  and  p(i7)  = 0,75/cos(i^), 

- Fig.  6d)  0,59  s & < 2,76  and  p{&)  - 0,5/sin(tf) . 

- Fig.  6e):  2,76  s < 4,02  and /y(i9)  = -l,25/cos(i9-) . 


Figure  6 Transform  of  the  strongly  singular  integrand 
into  regular  sub-domains  using  a transform  to 
polar  coordinates  and  a singularity  extraction. 
Integrands  in  sub-domains  are  normalized. 


3.3.2  Domain  transform  and  extraction  of 
singularity 


The  new  method  combines  the  domain  transform  and  the 
extraction  of  singularity  in  a convenient  way.  Analogous  to 
3.3.1,  the  strongly  singular  part  (9)  due  to  the  derivative  of 
the  Green's  function  is  treated.  The  order  of  the  singularity 
is  reduced  by  one  by  virtue  of  the  integration  weight 
duxdu2  = 2 ADtdsxds2  - AD‘(tx  + \)/4dtldt2  in  (19)  (see 

3.2,2).  This  weight  is  formed  by  the  transform  of  the 
differentials  from  the  w,,«2- system  into  the  tx  ,?2 -system. 

(15)  to  (18)  show  the  transform  for  wlt(f,,f2)  and  w2*(?i>0* 


k = 1...4  are  the  number  of  the  triangle  in  Fig.  4.  One 
obtains  for  (9)  with  (19)  for  the  x-component  the  weakly 
singular  integrand: 


*.’(*,  .0 


I 


An 


[-cos(jS'j?t(r,,f2))]A‘-^,t1)(38) 


In  (38),  the  singular  point  P is  mapped  onto  the  edge 
tx  ~ -1 . In  contrast  to  3.2.3  and  3.3.1 , the  distance  factor  is 
not  expressed  in  the  ux,u2 -system  around  the  point  P,  but 
rather  in  the  ft,f2 -system  along  the  edge  tx  = -1.  Using  a 
series  expansion,  the  extracted  term  is  analytically  integrable 
according  to  tx . For  the  x-component  (the  y-  and  z- 
components  are  analogous)  one  obtains: 


h+iMM2) 


(39) 

(40) 


with 


The  term  needed  for  the  extraction  of  the  singularity  is  given 
by  (39)  - (41).  With  (38),  they  yield: 
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1 y ad. 

1 dxk  ( tx , t2 ) 

16jt  H 

PK)  *, 

r,— 1 

1 V Ad‘ 

i 0 

16 n HA 

(',+!)//('>)  *. 

l 

The  terms  in  the  brackets  are  regular 


(42) 


- 4 l,-+l  (,-+l 

£ r ,/j 


i folks, 2) 


Xk{tj  ’^2  ) 


(43) 


dtldt2 


The  integration  is  carried  out  using  the  method  in  3.1 . 


For  the  second  term  in  (42),  the  Cauchy  limit  is  formed 
with  the  substitution  t = f,  + 1 : 


dt 


dtdt2 


For  tt  -►  0 , one  introduces  an  e -vicinity  as  in  (40) 

^ = e’  = '.7l0  + o('.3). 

one  then  obtains  from  (44) 

— t-lim  yA°‘ 

1 6tv  £_*°  4d 


1 dxk  (f,  ,t2  )| 


(44) 


(45) 


dt, 


(ln(2/t  (t2 ))  - In  ejdt2  (46) 


■ i_  VA°.'7  1 d<M.  in(2 fk(t2))dt. 

16*  £ J,WJ  1 1 ’’  ' 


with  (see  (37)) 

Urn  Y A0,  V 1 


df2  ln(f)  = 0 


(47) 


The  integration  method  is  now  visualized  by  the  example  of 
the  integrand  in  (43)  as  derived  from  (9).  The  strongly 
singular  integrand  in  the  wl,M2-system  is  in  Fig.  7a).  The 

integrands  after  transform  and  extraction  are  shown  in  the 
t,,t2-system  in  Fig.  7b)  to  7e).  The  strongly  singular  part 

in  (9)  is  transformed  into  a very  smoothly  varying  function. 
Hence,  the  integration  is  carried  out  accurately  with  only  a 
few  nodes  ensureing  a fast  computation. 


Figure  7 Transform  of  the  strongly  singular  integrand 
into  regular  sub-domains  using  a domain 
transform  and  a singularity  extraction. 
Integrands  in  sub-domains  are  normalized. 

3.4  Nearly  Singular  Integration 

The  integrands  in  Method  of  Moment  computations  are 
nearly  singular  when  structures  with  a small  thickness  are  to 
be  computed,  or  if  fields  very  close  to  a surface  must  be 
evaluated.  In  the  first  case,  the  problem  is  noticeable  when 
setting  up  the  matrix,  in  the  second  case  when  computing 
the  fields  after  the  surface  currents  have  been  computed. 
Figures  8a)  and  8b)  show  the  two  cases.  In  Fig.  8a),  the 
distance  between  two  disrcrete  surface  elements  is  very  small 
compared  to  the  dimensions  of  the  elements.  In  Fig.  8b), 
the  distance  d of  an  observation  point  P to  the  surface  is 
smaller  than  the  geometrical  dimensions  of  the  discretization 
elements. 


Figure  8 Cases  when  nearly  singular  integrals  are 
encountered. 


The  straight-forward  method  to  evaluate  nearly  singular 
integrals  is  to  increase  the  number  of  nodes  within  the 
numerical  quadrature  method.  In  this  case,  the  number  of 
nodes  increases  very  quickly  when  the  distance  d is  reduced. 

A second  method  is  to  solve  the  nearly  singular  integrals  by 
subdividing  the  surface  elements  in  smaller  elements,  and  to 
continue  the  subdivision  process  until  convergence  is 
achieved.  This  ensures  that  the  number  of  nodes  is  increased 
in  the  vicinity  of  the  observation  point.  The  subdivision 
process  is  numerically  expensive  and  is  inefficient  for  very 
small  distances  d.  Furthermore,  the  accuracy  of  a Gauss- 
Legendre  quadrature  depends  on  the  number  of  nodes:  The 
integration  using  5 nodes  of  an  element  that  has  been 
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subdivided  into  two  elements  allows  a correct  integration  of 
a polynomial  of  the  9th  order  on  each  subelement.  Without 
the  subdivision  process,  a polynomial  of  the  19th  order 
could  be  integrated. 

The  efficient  approach  in  [Tel87]  und  [Hay92]  is  to 
transform  the  surface  variables  ,w2  in  such  a way  that  the 

resulting  Jacobian  lifts  the  nearly-singularity.  This  is  to  say 
that  the  Jacobian  should  approach  the  distance  d at  the 
critical  point  [Tel87],  [Hay92].  The  transform  also 
concentrates  the  nodes  around  the  location  of  the  near- 
singularity which  leads  to  a better  integration. 

4.  Example  Integrals 


The  previously  presented  methods  for  integrations  in 
Method  of  Moments  computations  are  now  applied  to  an 
analytical  problem.  The  integrals  in  Table  3 are  derived  from 
the  integrals  in  Table  1 to  yield  the  analytically  integrable 
functionals  (48)  and  (49).  The  geometry  for  the  verification 
is  a quadratic  surface  element  that  is  centered  at  the  origin  as 
in  Fig.  1 . Each  side  has  a length  of  2 m at  a wavelength  of 
20  m.  Hence,  the  length  of  each  side  is  X/10.  The 
observation  point  P is  located  at  x = 0.25  m,  y - 0.5  m and 
a variable  height  z.  The  source  point  is  located  at  z = 0. 


Table  3 Example  integrals  for  the  verification  of  the 
integration  methods. 


1 rr  r'-rp 

weakly  singular  kernel  — JJ  rz - — 

4zr  Sr  |r  rp 

jdS'  (48) 

R -*  0 weakly  singular 

->  4.2 

1 rr  r>  _ rp 
strongly  singular  kernel  * — J J , r r „ 

4*  s'  r _ r 

|T dS'  (49) 

R > d 1 

— 4.1 

R < d J 

regular 

R -*0 

strongly  singular 

— 4.3 

R < d 

nearly  singular 

->  4.4 

The  integral  is  solved  analytically  and  numerically,  the 
following  error  evaluates  the  quality  of  the  method: 


error  in  % = 


[analytical  - numerical] 
|analytical| 


100 


(50) 


Note  that  it  is  only  possible  to  analytically  solve  the 
integrations  in  Table  3 for  planar  geometries,  not  for 
arbitrary  ones. 


Figure  9 Regular  integrand.  Error  for  various  locations  z 
of  observation  point  and  number  of  nodes. 

4.1  Regular  Integrals 

The  observation  point  P is  located  at  distances  from 
zP  = 0,01  m to  zP  - 2 m from  the  surface  element. 
Consequently,  the  kernel  varies  from  nearly  singular  R < d 
to  regular  R>  d.  An  integration  using  the  Gauss-Legendre 
in  3.1  with  nxn  nodes  yields  the  error  in  Fig.  9.  For  large 
distances  from  the  surface  to  the  observation  point  P,  very 
few  nodes  are  sufficient  for  a very  low  error.  The  error 
increases  with  decreasing  distance  z but  can  be  decreased  by 
using  more  nodes.  The  method  fails  for  very  small  distances 
z < 0.25  m. 

4.2  Weakly  Singular  Integrals 

The  observation  point  is  located  at  zP  = 0 on  the  surface  of 
the  element. 

4.2.1  Transform  to  Polar  Coordinates 

The  integration  of  (48)  is  carried  out  by  subdividing  the 
surface  element  into  four  sub-elements,  shown  in  Fig.  10 
(middle)  as  dashed  lines.  The  surface  integration  on  each 
quadrangle  composed  of  two  triangles  Dk  and  Di+1  with 

k=  1...8  (Fig.  10  rightmost)  again  involves  n x n nodes. 
til  2-nodes  are  used  for  the  tf-integration  per  triangle  Dk , the 

integration  according  p is  then  carried  out  using  n nodes. 


Figure  10  Surface  element  (z  = 0)  and  representation 
with  normalized,  local  coordinates.  The 
rectangle  is  subdivided  into  four  quadrangles 
composed  of  two  triangles  each. 
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The  integration  of  (48)  is  itemized  in  Table  4.  Convergence 
is  fast  with  increasing  number  of  nodes  to  the  correct  value. 
For  n > 32,  no  improvements  are  observed,  the  result  is 
stable  with  an  increasing  number  of  nodes. 

Table  4 Weakly  singular  integration  of  (48)  with 
transform  into  polar  coordinates.  Analytical 
result  is  -0.0579393378. 


n 

result 

error  (%) 

2 

-0.0519166274 

10.39 

4 

-0.0570911480 

14.64- 10-' 

6 

-0.0578289883 

19.05 -10-2 

8 

-0.0579257050 

23.53 -10'3 

16 

-0.0579393352 

45.73 -lO'7 

32 

-0.0579393378 

71.31-10-“ 

64 

-0.0579393378 

26.96 -10'10 

a) 


Figure  11  Two  discretizations  of  the  geometry  in  Fig.  1 
for  P(x  = 0.25  m,  y = 0.5  m,  z = 0). 


4.2.2  Domain  Transform 

The  surface  element  is  sub-divided  into  four  triangles  that 
are  in  turn  mapped  to  quadrangles.  Again,  the  integration 
involves  n x n-nodes  corresponding  to  the  number  of  nodes 
of  the  previous  example.  Table  5 shows  the  results  for  the 
method.  Again,  convergence  to  the  correct  value  is  observed 
with  a stable  result  beyond  n = 32  nodes. 


Figure  12  Surface  element  (z  = 0)  and  representation 
with  transform  to  polar  coordinates  (middle) 
and  domain  transform  (rightmost).  Polar 
coordinates:  The  rectangle  is  subdivided  into 
four  quadrangles.  Domain  transform:  The 
rectangle  is  subdivided  into  four  quadrangles 
composed  of  two  triangles  each.  Only  element 
1 (upper  right  comer)  is  shown. 


Table  5 Weakly  singular  integration  of  (48)  with 
domain  transform.  Analytical  result  is 
-0.0579393378. 


n 

result 

error  (%) 

2 

-0.0360523044 

37.78 

4 

-0.0549348239 

51.86  -10"1 

6 

-0.0578618391 

13.38  - 10-2 

8 

-0.0579760027 

63.28 -lO'3 

16 

-0.0579393184 

33.55-10-* 

32 

-0.5793933785 

22.88 -lO'10 

64 

-0.0579393378 

14.28-10-“ 

43  Strongly  Singular  Integration 

The  strongly  singular  integration  (49)  with  the  observation 
point  P(jc  = 0,25  m,  y = 0,5  m,  z = 0)  is  carried  out 
according  to  the  Cauchy  singular  integral.  The  integration 
domain  is  subdivided  into  four  sub-domains  as  in  Fig.  11. 
This  permits  the  evaluation  of  the  accuracy  of  the  surface 
integral  and  the  line  integral.  Fig  11a)  shows  a discretization 
into  only  one  surface  element.  Fig.  lib)  discretizes  the 
geometry  into  four  elements,  the  observation  point  is  then 
on  the  comer  common  to  each  discretization  element.  For 


the  discretization  in  Fig.  11a)  the  surface  integral  in  (49) 
vanishes  (see  3.3),  the  integration  is  completely  determined 
by  the  line  integral.  In  Fig.  lib)  both  terms  contribute  to 
the  result. 

The  transforms  in  turn  are  sketched  in  Fig.  12. 

4.3.1  Transform  to  Polar  Coordinates  and 
Extraction  of  Singularity 

[GG90]  evaluates  the  integral  for  many  combinations  of 
observation  points  and  surface  elements.  These  cases  have 
been  verified  with  the  present  method.  The  transform  is 
shown  in  Fig.  12  (middle). 

For  a discretization  into  one  element  (cf.  Fig.  11a),  the 
results  are  summarized  in  Table  6.  A good  convergence  is 
observed  for  an  increasing  number  of  nodes.  For  n = 6,  an 
error  of  almost  0.05  % is  achieved.  The  method  remains 
stable  for  a larger  number  of  nodes. 

The  integration  according  to  Fig.  lib)  maps  the  observation 
point  to  ux  - -I,  w2  « - 1 in  local  coordinates  (Fig.  12 

middle).  For  each  triangle,  n x n nodes  are  used,  the  line 
integral  uses  again  n nodes.  Table  7 shows  the  results  of  the 
integration.  More  nodes  are  needed  for  the  same  accuracy 
compared  to  a discretization  into  one  element  only. 
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Table  6 Strongly  singular  integration  of  (49)  after 
transform  to  polar  coordinates  according  to 
Fig.  10  and  11a).  Analytical  result  is 
-0.052741730991244. 


n 

surface  integral 

line  integral 

error  (%) 

2 

0 

-0.0471475986 

10.61 

4 

0 

-0.0524191674 

61.16  10"2 

6 

0 

-0.0527116953 

56.95  10-3 

8 

0 

-0.0527389852 

52.06-10^ 

16 

0 

-0.0527417307 

49.36  10-8 

32 

0 

-0.0527417310 

54.90- 10-“ 

64 

0 

-0.0527417310 

18.89*10'° 

Table  7 Strongly  singular  integration  of  (49)  after 
transform  to  polar  coordinates  according  to 
Fig.  10  and  12.  Analytical  result  is 
-0.052741730991244. 


n 

surface  integral 

line  integral 

error  (%) 

2 

0.2416862365 

-0.1793554832 

21.82  *10+1 

4 

0.1785239348 

-0.2481406072 

32 

6 

0.1849263635 

-0.2398825203 

41.99  -10"1 

8 

0.1868837023 

-0.2388539817 

14.63*  10'1 

16 

0.1864218113 

-0.2391635614 

36.29-10'6 

32 

0.1864218178 

-0.2391635488 

12.23  *10-9 

64 

0.1864218178 

-0.2391635488 

19.75  *10'10 

Table  8 Strongly  singular  integration  of  (49)  by 
domain  transform  according  to  Fig.  11a). 
Analytical  result  is  -0.052741730991244. 


n 

surface  integral 

line  integral 

error  (%) 

2 

0 

-0.0883362945 

67.49 

4 

0 

-0.0615749331 

16.75 

6 

0 

-0.0529196853 

33.74  -10'2 

8 

0 

-0.0525051309 

44.86-1  O'2 

16 

0 

-0.0527419384 

39.33  - 10-5 

32 

0 

-0.0527417310 

91.62-  lO-'1 

64 

0 

-0.0527417310 

22.17 -10-10 

Table 

9 Strongly  singular  integration  of  (49)  b> 
domain  transform  according  to  Fig.  12 
Analytical  result  is  -0.052741730991244. 

n 

surface  integral  line  integral 

error  (%) 

2 

0.2330166165  -0.1652108191 

22.86-10*' 

4 

0.18801 15487  -0.2448561289 

77.79  -lO'1 

6 

0.1809303725  -0.242759201 1 

17.23 

8 

0.1869992436  -0.2386698415 

20.31-10-' 

16 

0.18642321 12  -0.2391626624 

43.23 -10-4 

32 

0 . 1 8642 1 8 1 78  -0 .239 1 63548  8 

19.17-10"9 

64 

0.1864218178  -0.2391635488 

39.71-10-'° 

4.3.2  Domain  Transform  and  Extraction  of 
Singularity 

The  integration  with  domain  transform  and  subsequent 
extraction  of  the  singularity  yields  the  results  in  Table  8. 
The  integration  shows  stable  results  for  an  increasing 
number  of  nodes.  An  accuracy  of  better  than  0.01  % is 
reached  for  16  nodes. 

A discretization  according  to  Fig.  lib)  results  in  a sub- 
division as  shown  in  Fig.  12  for  the  element  1.  The 
singularity  is  located  at  w,  = -1,  u2  = -1  for  each  triangle. 

The  integration  is  carried  out  according  to  3.3.2,  n x n 
nodes  are  used.  Table  9 shows  the  results  for  both  surface 
and  line  integrals.  A stable  convergence  is  achieved,  an  error 
better  than  0.01  % is  obtained  for  n = 16. 

The  integrations  in  3.3.1  and  3.3.2  yield  both  results  within 
the  same  range  of  accuracy.  For  the  result  shown,  the 
method  in  3.3.2  (viz.  the  domain  transform  and  subsequent 
extraction  of  singularity)  results  in  smaller  errors  for  only 
few  nodes. 


4.4  Nearly  singular  integrals 

As  seen  in  3.1,  the  integration  using  a Gauss-Legendre 
quadrature  when  the  observation  point  is  close  to  a 
surface  (see  Fig.  1)  yields  only  slowly  converging  results. 
A large  number  of  nodes  is  required.  Even  for  n = 32,  the 
error  reaches  3000  % for  an  observation  point  located  at 
zP  =0.01  m . The  method  in  [Tel 87],  [Hay 92]  allows  a fast 
integration  with  a high  level  of  accuracy. 

Fig.  13  shows  the  error  obtained  with  the  method  in  3.4  for 
an  increasing  number  of  nodes.  The  shortest  distance 
between  an  observation  and  a source  point  is  zP  = 0.01  m . 
The  largest  error  is  for  n = 16  at  a distance  of  zP  - 0.02  m 
and  evaluates  to  2.1631  %.  For  n = 32,  this  error  is  reduced 
to  0.2106  %. 

Figure  14  shows  a comparison  between  the  methods  in  3.1 
and  3.4  for  the  same  number  n = 4 and  n = 32.  It  is  clear 
from  this  picture  that  the  integration  according  3.4  yields 
much  better  results  for  nearly  singular  integrals. 
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Figure  13  Nearly  singular  integrand.  Error  curves  for  an 
increasing  number  of  nodes  with  observation 
point  at  various  distances  to  the  surface 
element. 

5*  Conclusions 

By  properly  identifying  the  orders  of  the  singularity  with  the 
reaction  integrals  of  the  Method  of  Moments  formulation, 
appropriate  integration  methods  are  found  to  be  applicable 
for  each  of  the  integrals.  Regular  integrals  are  solved  by  a 
Gauss-Legendre  quadrature.  Weakly  singular  integrals 
necessitate  a regularization,  and  are  then  integrated 
numerically.  Strongly  singular  integrals  are  regularized  by  a 
two  stage  process  and  integrated  both  analytically  and 
numerically.  The  correct  method  ensures  for  each  kernel  a 
fast  and  stable  convergence  to  the  correct  values. 

The  methods  are  general  enough  to  be  used  as  methods  for 
surface  integrals  with  various  degrees  of  singularity.  The 
methods  are  not  restricted  to  Method  of  Moments 
formulations  in  Computational  Electromagnetics. 
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